A new automatic geo-electric self-potential imaging technique for diverse sustainable development scenarios

This study introduces a rapid and efficient inversion algorithm designed for the interpretation of self-potential responses originating from mineralized and ore sources and hydrothermal activity, specifically addressing spherical, vertical, and horizontal cylindrical structures. The algorithm leverages local wavenumber and correlation imaging techniques to enhance accuracy in modeling. The correlation factor (Cf value) is crucial in this approach, calculated as the correlation between the local wavenumber of the measured self-potential field and that of the computed field. The algorithm identifies the maximum correlation Cf value (CF-max) as indicative of the optimal true model parameters. To validate the proposed algorithm, it was applied to three theoretical examples—one with contamination from regional background and another with multiple sources with and without different types of noises (random Gaussian and white Gaussian noises). Additionally, the approach was tested on three distinct real field cases related to mining, ore investigation and hydrothermal activity in India, Germany and USA. Through a comprehensive analysis of results from theoretical and real-world scenarios, including comparisons with different available data and literature information, the study concludes that the method is effective, applicable to multiple sources, accurate, and does not necessitate prior knowledge of the source shape. This algorithm presents a promising advancement in the field of self-potential interpretation for mineral exploration and geothermal exploration.


Methodology
The self-potential signature (P) at an observation point (x j , z) along the profile depicted in Fig. 1, can be expressed using the formula of Yüngül 39 .
where n represents the count of data points, q is the shape factor, a dimensionless quantity, varies according to the structure's shape (it takes a value of 1.5 for a structure resembling a sphere, 1 for a horizontal cylindrical structure, and 0.5 for a semi-infinite vertical cylinder.The depth of the structure is denoted as 'z o ' in meters.The amplitude factor (K), with unit mVm 2q−1 , the parameter 'x o ' indicates the position of the source body in meters, and ' θ ' corresponds to the polarization angle in degrees.
The measured local wavenumber can be formulated by 30,40 : where (1) P K, x j , x o , z o , θ, q = K (x l − x o )cosθ + z o sinθ www.nature.com/scientificreports/By substitution (Eq. 3 in Eq. 2), LW mea can be given by: where AS is the analytical signal amplitude as follow 41 : By applying the horizontal and vertical derivatives ∂P ∂x and ∂P ∂z respectively to Eq. (1) and substituting in Eq. ( 4), the computed local wavenumber ( LW com ) is given by: Using LW mea and LW com , the correlation parameter can be represented by 30,40 : Using Eq. ( 7), the calculation of the correlation parameter ( C F ) between LW mea and LW com is performed, and the highest value of C F corresponds to actual body characteristics 30,40 .The process flow of the proposed algorithm is illustrated in Fig. 2.After identifying the most suitable parameters from the search space based on the highest C F value (C F -max), it becomes possible to create a two-dimensional representation of C F for the preferred source (specifically, the shape factor q) in relation to subsurface depth (m).The solid black dot present in the imaging section symbolizes the accurate position for both depth and location.

Synthetic models
This segment demonstrates the application of the suggested method to three distinct synthetic models, both with and without noise, in order to assess the effectiveness and suitability of the proposed approach in the interpretation of self-potential anomalies.

Example 1
The self-potential profile resulting from a horizontal cylinder was computed with specific parameters: K = 3500 mV m, z o = 10 m, x o = 0 m, q = 1, and θ = -55°, over a profile length of 100 m (depicted in Fig. 3a).The interpretation process began by calculating both horizontal and vertical gradients of the observed anomaly (as shown in Fig. 3b).Subsequently, the value of LW mea was determined using Eq.(4) (illustrated in Fig. 3c).Moving forward, the calculation of C f was carried out using Eq. ( 7) (as demonstrated in Fig. 3d), considering various q values as presented in Table 1.Notably, in Table 2, the highest value of C f (C F -max = 1) (black circle in Fig. 3d) is located at K = 3500 mV m, z o = 10 m, x o = 0 m, q = 1, and θ = -55°, which aligns with the information in Fig. 3d.This outcome signifies the exceptional efficiency of the proposed method.Utilizing the suggested approach facilitated the estimation of inverted parameters as detailed in Table 2, leading to a complete absence of errors for the diverse parameters.
Secondly, 15% WGN (Fig. 5a), the noisy data's vertical and horizontal gradients were computed (Fig. 5b).Subsequently, applying Eq. ( 4), LW mea was determined (Fig. 5c).To derive C f , Eq. ( 7) was employed (Fig. 5d).Within Fig. 5d, the highest C f value of 0.6784 (depicted by the black circle in Fig. 5d) is observed at K = 4015 mV m, z o = 12 m, x o = -− 1 m, q = 1, and θ = − 54.5°, as indicated in Table 3.The computed error of the estimated parameters, K, z o , q, θ are: 14.7%, 20%, 0% and 0.91% respectively.The results obtained above shows that the effect of the WGN is greater than RGN on the proposed method but the estimated parameters in case of the different two types of noise demonstrating that the proposed method can effectively be employed to handle noisy data with exceptional performance. (3)

Example 2
In order to evaluate the suitability and effectiveness of the employed approach when dealing with multisource examples, the technique was implemented on a 100 m composite profile that was constructed of sphere body (applying these parameters: K = 30,500 mV m 2 , z o = 5 m, x o = − 30 m, q = 1.5, and θ = − 25°) and horizontalcylinder (HC) body (applying these specific parameters: K = 2500 mV m, z o = 3 m, x o = 30 m, q = 1, and θ = − 35°) (Fig. 6a).The process of interpretation began with the computation of both the vertical and horizontal gradients of the composite profile (depicted in Fig. 6b).The value LW mea was determined applying Eq. (4) (illustrated    Table 3.The authentic and retrieved model parameters pertaining to the first theoretical example (selfpotential anomaly generated by a horizontally cylinder) contaminated with 15% RGN and 15% WGN. in Fig. 6c), while the value of C f was calculated using Eq. ( 7) (as shown in Fig. 6d).It is demonstrated that the highest C f value (C F -max 1 ) for the sphere body is 0.74 (indicated by the initial black circle in Fig. 6d), situated at K = 27,718 mV m 2 , z o = 4.9 m, x o = − 30 m, q = 1.5, and θ = − 24.71° (Table 4).Similarly, the maximum C f value (C F -max 2 ) for the HC body is 0.70 (depicted by the second black circle in Fig. 6d), located at K = 2636 mV m, z o = 3.4 m, x o = 30 m, q = 1, and θ = − 38.79°, as summarized in Table 4.The computed error of the estimated parameters, K, z o , x o , q, and θ are: 9.12%, 2%, 0%, 0% and 1.16% respectively for the sphere body, while for the HC source, the computed error of the estimated parameters, K, z o , x o , q, and θ are: 5.44%, 13.33%, 0%, 0% and 10.83%, respectively.

Recovered
To evaluate how well the proposed method performs in the presence of noise and its overall effectiveness, the method was applied to the previous model after adding 10% RGN and 10% WGN.The first case, 10% RGN (Fig. 7a), the noisy composite data's vertical and horizontal gradients were computed (Fig. 7b), The value LW mea was determined (illustrated in Fig. 7c), while the value of C f was calculated (as shown in Fig. 7d).It is demonstrated that the highest C f value (C F -max 1 ) for the sphere body is 0.55 (indicated by the initial black circle in Fig. 7d), situated at K = 23,287 mV m 2 , z o = 4.6 m, x o = − 30 m, q = 1.5, and θ = − 26.01° (Table 5).Similarly, the maximum C f value (C F -max 2 ) for the HC body is 0.56 (depicted by the second black circle in Fig. 7d), located at K = 2293 mV m, z o = 3.2 m, x o = 30 m, q = 1, and θ = − 42.67°, as summarized in Table 5.The computed error of  the estimated parameters, K, z o , x o , q, and θ are: 23.65%, 8%, 0%, 0% and 4.04% respectively for the sphere body, while for the HC source, the computed error of the estimated parameters, K, z o , x o , q, and θ are: 8.28%, 6.67%, 0%, 0% and 21.91%, respectively.The second case, 10% WGN (Fig. 8a), the noisy composite data's vertical and horizontal gradients were computed (Fig. 8b), The value LW mea was determined (illustrated in Fig. 8c), while the value of C f was calculated (as shown in Fig. 8d).It is demonstrated that the highest C f value (C F -max 1 ) for the sphere body is 0.50 (indicated by the initial black circle in Fig. 8d), situated at K = 36,322 mV m 2 , z o = 5.9 m, x o = − 30 m, q = 1.5, and θ = − 27.55° (Table 5).Similarly, the maximum C f value (C F -max 2 ) for the HC body is 0.48 (depicted by the second black circle in Fig. 8d), located at K = 2833 mV m, z o = 4.3 m, x o = 30 m, q = 1, and θ = − 47.50°, as summarized in Table 5.The computed error of the estimated parameters, K, z o , x o , q, and θ are: 19.09%, 18%, 0%, 0% and 10.2% respectively for the sphere body, while for the HC source, the computed error of the estimated parameters, K, z o , x o , q, and θ are: 13.32%, 43.33%, 0%, 0% and 35.71%, respectively.
Hence, it can be inferred that the suggested approach is well-suited for scenarios involving multiple sources.

Example 3
To assess the effectiveness of our approach when dealing with a regional context, we introduced a self-potential anomaly profile originating from a vertical cylinder (characterized by: K = 250 mV, z o = 4 m, x o = − 25 m, q = 0.5, and θ = − 75 o and profile length 100 m) into a deep-seated first order regional anomaly (Fig. 9a).The interpretation process commenced by computing both the horizontal and vertical gradients of the observed anomaly, as depicted in Fig. 9b.Subsequently, Eq. ( 4), was employed to ascertain the value of LW mea (Fig. 9c).
To evaluate the robustness and efficacy of the proposed method when applied to noisy data, we introduced two types of noise, specifically, 15% RGN and 15% WGN, to the previous model.

Field models
In order to evaluate the effectiveness of the suggested method, it was employed in three distinct real-life field data, including one from India, one from Germany and the third from USA.

India field example (Neem-Ka-Thana Copper Belt)
The Neem-Ka-Thana Copper Belt in India is distinguished by the prevalence of copper mineralization in the region 4,42 .Notably, the copper deposits are primarily located along fault lines and shear planes, indicating a geological association with these structural features.The concentration of copper in the mines within the Neem-Ka-Thana Copper Belt shows variability, ranging from 0.6 to 1.2% 4,43 .This diversity in copper content underscores the geological complexity of the region, suggesting that mineralization processes have been influenced by a combination of tectonic forces and geological phenomena.Therefore, the Neem-Ka-Thana Copper Belt stands out as a significant geological site where the interplay of geological structures and mineralization processes contributes to the formation of valuable copper deposits.A profile of self-potential was acquired over the Neem-Ka-Thana Copper Belt in India 4,31,43 (Fig. 12a).The profile length was 285 m long.The interpretation procedure initiated with the computation of both the horizontal and vertical gradients of the observed anomaly, illustrated in Fig. 12b.Following this, Eq. ( 4) was applied to determine the value of LW mea , as shown in Fig. 12c.Progressing further, Eq. ( 7) was employed to compute C f , depicted in Fig. 12d, considering various q values as presented in Table 8.It is noteworthy that Table 8 presents the maximum value of C f (C F -max = 0.96), represented by a black circle in Fig. 12d, corresponding to K = − 47.92 mV, z o = 18 m, x o = 177.5 m, q = 0.4, and θ = 88° (Table 9).The comparison between the results obtained by our suggested method and those obtained by Balkaya 44 is depicted in Fig. 12a.Table 9 displays a comparison of the inverted parameters between the proposed method and those of different methods found in the literature.

Germany field example (Lias-epsilon black shales)
The Lias-epsilon black shales in Germany are situated atop a coal maturity high on the Bramsche Massif in Northwest Germany, as described by 45 .The thermal evolution of this region is attributed to the inversion of the Lower Saxony Basin, occurring during the Early Late Cretaceous period, likely in conjunction with mafic intrusions from the Bramsche, Vlotho, and Uchte Massifs at depths of approximately 5-10 km [46][47][48] .Table 6.The authentic and retrieved model parameters pertaining to the second theoretical example (selfpotential composite anomaly induced by vertical cylinder and first order regional source).www.nature.com/scientificreports/ The heightened thermal maturity of organic materials in the Lias-epsilon black shales is commonly associated with the intrusion of the Bramsche Massif into the Earth's crust.Notably, strong gravity and magnetic anomalies, along with increased coal maturity in the Westphal D coals found in areas such as Ibbenbüren's mining region, are linked to the presence of the Bramsche Massif [49][50][51][52] .The thermal heating of the stratigraphic series likely commenced before the Alp era and extends beyond Mesozoic black shales (2-5% C-org) of the Lower Toarcium and Lias-epsilon, as indicated by Mann 53 .The contemporary morphology of the region has been significantly influenced by Pliocene tectonism and Quaternary sedimentation, as highlighted by studies such as those conducted by 54,55 .
The survey area's location is depicted in Fig. 13 45,49 .A self-potential profile was carried out across a 500 m span over the conductivity anomaly, specifically the Lias-epsilon black shales 45,56 (Fig. 14a).The interpretation process began by calculating both the horizontal and vertical gradients of the observed anomaly, as illustrated in Fig. 14b.Subsequently, Eq. ( 4) was utilized to determine the value of LW mea , as depicted in Fig. 14c.Advancing further, Eq. ( 7) was applied to calculate C f , shown in Fig. 14d, with consideration for various q values outlined in Table 10.It is important to note that Table 10 highlights the maximum value of C f (C F -max = 0.71), denoted by a black circle in Fig. 14d.This corresponds to K = 11,052.05mV m, z o = 19 m, x o = 250 m, q = 1, and θ = − 100° (refer to Table 11).The comparison between the results obtained by our suggested approach and those obtained   57 is depicted in Fig. 14a.Table 11 provides a comparison of the inverted parameters between the proposed method and those from different methods documented in the literature.Also, Figs. 15 and 16 show the depth estimation from the 2D electrical resistivity tomography and the estimated source model using the suggested technique, respectively (taking into consideration the topography of the area).The results from the 2D electrical resistivity tomography and our method match well.

USA field example (Hi'iaka eruption)
On May 5, 1973, a dike penetrated the upper crust of Kilauea volcano within the geologic context of the east rift zone 58 .This intrusion coincided with the eruption of Hi'Iaka and Pauahi craters, as documented by Klein et al. 59 and Tilling et al. 60 .The dike induced the formation of a surface fissure, stretching 100 m, which erupted magma west-southwest (WSW) of Hi'iaka crater.Geophysical measurements indicated that the dike extended underground in the WSW direction for an additional 1.5 km 58 .
The SP profile's location is depicted in Fig. 17 58 .The selected profile was carried out in 1997 58 , spanning a length of 650 m (see Fig. 18a).The process of interpretation commenced by computing both the horizontal and vertical gradients of the observed anomaly, as depicted in Fig. 18b.Subsequently, Eq. ( 4) was applied to ascertain the value of LW mea , as illustrated in Fig. 18c.Progressing further, Eq. ( 7) was employed to compute C f , as shown in Fig. 18d, considering various q values outlined in Table 12.It is noteworthy that Table 12 highlights the maximum value of C f (C F -max = 0.89), represented by a black circle in Fig. 17d.This corresponds to K = − 4688.45 mV m 2q-1 , z o = 110 m, x o = 320 m, q = 0.7, and θ = − 110° (refer to Table 13).The comparison between the results obtained by our suggested method and those obtained by Mehane et al. 57 is depicted in Table 9. Retrieved model parameters for the Neem-Ka-Thana Copper Belt field example in India, and the comparison of the inverted results between the proposed method and those of different methods found in the literature.

Model parameters
Balkaya 44 Göktürkler and Balkaya 31 Biswas 4 13 presents a comparison of the inverted parameters between the proposed method and those documented in the literature from different methods.

Conclusions
In this study, we implemented an effective inversion imaging algorithm to characterize self-potential data originating from diverse sources such as spheres, vertical cylinders, and horizontal cylinders.The demonstrated algorithm holds promise for applications in mineral, ore exploration, and geothermal investigation offering precise predictions of various structural parameters-namely, amplitude factor (K), depth (z o ), body origin (x o ), shape factor (q), and polarization angle ( θ)-with high accuracy and without the need for a priori information.The sug- gested algorithm employs the correlation factor (C f ) between the local wavenumber of the observed self-potential field and that of the computed field.The findings indicate that the maximum C f (C F -max) corresponds to the most reliable estimated model.Moreover, our proposed approach presents an imaging algorithm that provides rapid (within seconds) and robust imaging for subsurface depth and the location of concealed anomalous sources.To validate the efficiency, accuracy, and stability of the proposed algorithm, we subjected it to testing using three synthetic cases, including a pure data, a noisy data contaminated with different types of noise (RGN and WGN), an example for multi-source model and data with regional background effects.The applicability of the algorithm was further assessed through three real cases for mineral/ore exploration and geothermal investigation in India, Germany and USA.The resulting models from these real cases exhibited strong correlations with drilling data and findings reported in the literature.Finally, our study supports the suitability of the proposed algorithm for mineral/ore deposits exploration and geothermal investigation as well.
Table 11.Retrieved model parameters for Lias-epsilon black shales field example in Germany, and the comparison of the inverted results between the proposed method and those of different methods found in the literature.
Model parameters Gurk et al. 45 Mehanee 61 Mehanee et al. 57 Present study K (mV m       Table 13.Retrieved model parameters for the Hi'iaka eruption field example in USA, and the comparison of the inverted results between the proposed method and those of different methods found in the literature.
Model parameters Davis 58 Mehanee et al. 57 Present study K (mV m 2q

Figure 2 .Figure 3 .
Figure 2. Flowchart depicting the procedural sequence of the algorithm under consideration.

Figure 4 .
Figure 4. (a) Profile of the self-potential anomaly depicted in Fig. 3a after contaminating with 15% RGN, (b) computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (d) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

15%Figure 5 .
Figure 5. (a) Profile of the self-potential anomaly depicted in Fig. 3a after contaminating with 15% WGN, (b) computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (d) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

Figure 6 .
Figure 6.(a) Profile of the self-potential anomaly induced by multisource models of sphere and horizontal cylinder, (b) computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (c) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

Table 4 . 2 zFigure 7 .
Figure 7. (a) Profile of the self-potential anomaly depicted in Fig. 6a after contaminating with 10% RGN, (b) computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (d) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

Figure 8 .
Figure 8.(a) Profile of the self-potential anomaly depicted in Fig. 6a after contaminating with 10% WGN, (b) computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (d) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

Figure 9 .
Figure 9. (a) Profile of the self-potential composite anomaly induced by vertical cylinder and first order regional source, (b) computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (d) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

3 Figure 10 .
Figure 10.(a) Profile of the self-potential anomaly depicted in Fig. 9a after contaminating with 15% RGN, (b) computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (d) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

Figure 11 .Figure 12 .
Figure 11.(a) Profile of the self-potential anomaly depicted in Fig. 9a after contaminating with 15% WGN, (b) computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (d) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

Figure 13 .Figure 14 .
Figure 13.Location map of the of the survey area for the Osnabrück in Germany (after Stadler and Teichmuller 49 and Gurk et al. 45 ).NL Netherlands, B Belgium.
Fig. 18a.Table13presents a comparison of the inverted parameters between the proposed method and those documented in the literature from different methods.

Figure 15 .
Figure 15.Results of the 2D electrical resistivity tomography inversion cross-section for the Osnabrück anomaly in Germany.(modified from Gurk et al. 45 ).

Figure 16 .
Figure 16.Estimated subsurface model using our suggested technique for the Osnabrück anomaly in Germany, taking into consideration the surface topography (b).The plus sign indicates the center location of source anomaly (a).

Figure 17 .
Figure 17.Location map of the of the SP profile for the Hi'iaka eruption in USA (after Davis 58 ).

Figure 18 .
Figure 18.(a) Profile of the self-potential anomaly over Hi'iaka eruption in USA, along with the calculated responses from the present study and those obtained by Mehanee et al. 57 .(b) Computed horizontal and vertical derivatives for the profile depicted in (a), (c) local wavenumber of the data depicted in (b), (d) visualizing the correlation factor (C f ) and determining the C F -max through the newly established method.

Table 1 .
The correlation factor (C F ) calculated at the different shape factors for the first theoretical example (self-potential anomaly induced by a horizontal cylinder).The optimum values are in[bold].

Table 2 .
The authentic and retrieved model parameters pertaining to the first theoretical example (selfpotential anomaly generated by a horizontally cylinder).

Table 5 .
The authentic and retrieved model parameters pertaining to the second theoretical example (induced by multisource models of sphere and horizontal cylinder) contaminated with 10% RGN and 10% WGN.

Table 8 .
The correlation factor (C f ) calculated at the different shape factors for the Neem-Ka-Thana Copper Belt field example in India.The optimum values are in [bold].

Table 10 .
The correlation factor (C f ) calculated at the different shape factors for the Lias-epsilon black shales field example in Germany.The Optimum values are in [bold].

Table 12 .
The correlation factor (C f ) calculated at the different shape factors for the Hi'iaka eruption field example in USA.The Optimum values are in[bold].